Equation of motion method for Full Counting Statistics: Steady state superradiance 
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For the multi-mode Dicke model in a transport setting that exhibits collective boson transmissions, 
we construct the equation of motion for the cumulant-generating function. Approximating the 
exact system of equations at the level of cumulant-generating function and system operators at 
lowest order, allows us to recover master equation results of the Full Counting Statistics for certain 
parameter regimes at very low cost of computation. The thermodynamic limit, that is not accessible 
with the master equation approach, can be derived analytically for different approximations. 
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I. INTRODUCTION 

The fluctuations of particles tunneling randomly 
through a quantum structure are sensitive indicators for 
particle interactions [1], coherences [2] or collective ef- 
fects [3] . Originating from the field of quantum optics [4] , 
the theory of Full Counting Statistics (FCS) has proven 
to be a versatile tool for the evaluation of these fluctu- 
ations, that can be conveniently described by the cumu- 
lants of the underlying stochastic process [5]. FCS is 
successfully applied in electronic systems [6, 7], where 
the charge fluctuations can be measured by a nearby 
quantum point contact [8] . While for various small-sized 
quantum systems the FCS can be derived analytically [9] , 
the calculation of the full statistics, i.e., all higher cu- 
mulants, of the emitted particles for an arbitrarily large 
system size is cumbersome due to the required diagonal- 
ization of the large Liouvillian C(x), that describes the 
system dynamics. 

Apart from numerical methods the most obvious ap- 
proach to this problem is to make further approximations 
concerning the factorization of higher order correlation 
functions, as is routinely done in, e.g., the Hartrcc-Fock 
approach in solid state physics [10]. While this factor- 
ization is usually applied at the level of specific observ- 
ables, e.g., for two-electron processes, in this work we pro- 
pose a factorization for the correlations of the cumulant- 
generating function (CGF) and system operators derived 
within a weak-coupling theory. One of the most basic and 
yet size-scalable models is the Dicke model [11] for super- 
radiant decay of an initially excited atomic cloud inter- 
acting with a radiation field [12]. In the limit of small 
sample size, the dynamics of the system can be analyzed 
in the symmetrical angular momentum basis [13] and an- 
alytic results for higher correlations of the transient pulse 
of emitted bosons can be obtained [14]. A semi-classical 
propagator for transient Dicke superradiance was derived 
in [15]. 

In this work, we consider an extended multi-mode 
Dicke model in a transport setup [16], in a limit where 
back-action between counted particles and the system 
state can be neglected. Contrary to previous work we 
have therefore a steady-state transport setting that al- 



lows us to obtain analytic long-time results. For this 
aim, we construct an equation of motion (EoM) for the 
CGF, the solution of which requires us to make factoriza- 
tion assumptions for correlations of counting and system 
operators. We show that in certain limits even factoriza- 
tion at lowest level allows one to retrieve the FCS for ar- 
bitrary system size and access the thermodynamic limit. 
The semi-classical results for transient effects can be ob- 
tained as limiting cases [17]. 

With the recent progress in single-photon detectors [18] 
and the rising interest in collective boson transport, e.g., 
in thermal transport [19] or new types of lasers, like 
phonon and superradiant lasers [20, 21], we expect our 
findings to be relevant for a broader community. 

This work is structured as follows. In Sec. II we intro- 
duce the model under consideration and briefly introduce 
the underlying theory. In Sec. Ill we present exact EoM 
for the CGF and correlations. In Sec. IV we consider 
arbitrary system size N and compare the solutions for 
different approximations that close the EoM (Sec. IV A). 
We then derive analytic expressions for the thermody- 
namic limit of the FCS (Sec. IV B), examine the scaling 
behavior of the cumulants (Sec. IV C), and discuss the 
quality of the different approximations (Sec. IV D). We 
finish in Sec. V with a conclusion. 



II. MODEL 

Our model consists of two bosonic reservoirs, source S 
and drain D, that are coupled collectively by a medium 
of N two- level systems with identical level splitting f2. It 
is described by an extended multi-mode Dicke Hamilto- 
nian [16] 
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with collective spin operators J x z = X^i^f ' 2 ! where 
the operator b\ a creates a boson of frequency ujk a in 
reservoir a with the corresponding coupling constant hka- 
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FIG. 1. Sketch of the model: Bosonic excitations are ab- 
sorbed/emitted with rates Ts/d m and out of the medium of 
iV two-level systems to source or drain, respectively (occupa- 
tions ns > tid)- The number of total boson emissions into 
the drain is counted to construct the distribution P n (t), de- 
noting the probability for n particles effectively emitted into 
the drain after time t. 



We obtain the counting statistics of the total num- 
ber n of bosons exchanged between the medium and the 
drain by formally introducing [22] a bookkeeping oper- 
ator eft — X)^L-oo \ n 1) ( n l m th° Hamiltonian via 
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eft and similarly for the annihilation op- 
erator. The bookkeeping operator eft increases the occu- 
pation of a virtual detector by one unit for every boson 
created in the drain, independent on the mode k. In the 
following, we consider the TV-two-level medium and the 
bookkeeping operator as the system and assume a weak 
coupling between system and reservoirs. This allows us 
to derive a Lindblad master equation in Born, Markov, 
and secular approximation, where we assume initial de- 
coupling of system and bath, a memory-less bath and 
neglect fast rotating terms in comparison with the sys- 
tem time-scale [23]. 

With the extension of the Hamiltonian by the book- 
keeping operator, the system state = (n\ p \n) is 
conditioned on the number of bosons n measured in the 
detector. The master equation assumes the form 



(2) 



where C± describe emissions out of (+) and into (-) the 
system to and from the drain, while Co describes the 
remaining evolution. Explicitly, the super-operators are 
given by [16] 
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+T S [1 + ns] J-P-J+ + T s n s J+pJ- , 
C + p = T D [1 + n D ] J-pJ+ , 

C-p = T D n D J + pJ- , (3) 

where r a = 2irJ2k \hka\ 2 S(fl — iOka) denotes the sponta- 
neous boson emission rate of a single two-level system 
into a vacuum reservoir a and { • , • } is the anticommuta- 
tor. The numbers n a denote the stationary occupation 
of bath modes at system transition frequency ft: For the 



case of thermal baths one has, e.g. 
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The probability of counting n bosons after time t is 
given by P n (t) = Tr { p( n > (t) } . The n-resolved master- 
equation (2) can be Fourier-transformed by introducing 
a counting field \ Yia p(x>t) = Y^ n P ^ n \t)e mx leading 
to a generalized master equation for collective sponta- 
neous emission and absorption ^p(x,i) = C(x)p{Xit) 
with £(x) = £o + e +lx C+ + e~ IX £_. In case of a single 
reservoir and neglect of the counting field, this reproduces 
previous master equations [25]. 

The FCS of the bosonic detection probability P n (t) for 
n bosons after time t counted at the drain reservoir is 
given by the cumulants = ( — i<9 x ) fc C(x, £)| x =o, 

denoted by double brackets, where 

C( X ,t) = lnTr{p(x,*)} = InTr {e^po} (4) 

is the CGF, cf. e.g., [26]. In the stationary limit 
the CGF can be approximated by the eigenvalue Ao(x) 
with smallest modulus, fulfilling Ao(0) = 0, such that 
lim C(x,i) ~ Ao(x)^- Throughout this work we are in- 

t— > oo 

vestigating the long-time current cumulants ((Ik)) = 

lim f ^oo aj (("&(*)))• 

Calculations are conveniently performed in the angular 
momentum basis, where j = ^ and 

Jz | j, m) = 2m | j, m) , 



J± \j, m ) = VjO' + 1) - m(m±l) | j, m ± 1) , 
J 2 \j,m) =4j(j + l)\j,m) , (5) 

with J = (J x , J y , J Z ) T , and [J+, J_] = J z . In our model 
the system size N is fixed, so that we can treat the ex- 
pectation value / J 2 ) = J 2 = 2N(y + 1) as a constant. 
For later use note that J_ J + = | [j 2 — J 2 ] — |j z and 
[J„J±]=±2J±. 

III. EXACT RESULTS 

For the stationary properties of the expectation val- 
ues (J") we have an exact benchmark: Even for the 
case of non-equilibrated reservoirs the stationary state 
of the system, defined by C(0)p = 0, is a thermal one, 

p oc e~~ ~, with a temperature j3~ that corresponds to 
coupling the system to a single fictitious reservoir with 
the weighted average occupation 



T s n s + Tpn D 

r s + r D 
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where j3 a is the inverse temperature of bath a. 



Such simple average occupations exist also for an ar- 
bitrary number of baths whenever the operator struc- 
tures of the couplings are identical [24]. Since in addi- 
tion our system only supports a single allowed transition 
frequency, it effectively thermalizes at some average tem- 
perature even when coupled to baths at different temper- 
atures. 

Using this effective thermalization we find the ratio of 
the matrix elements p m = (jm\ p \jm) of the tri-diagonal 
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Liouvillian Pm+1 

Pm 



n+l ' 



In connection with normal- 



ization ^2 m Pm = 1) this implies that the steady state 
expectation values 



(7) 



with k = 1, 2, . . . can be calculated exactly (not shown 
for brevity) for arbitrary system sizes N. 

The evolution of the CGF for the emitted 
bosons at the drain can be calculated using 
C(x,t) = \n(e inx ) = \nJ2 n e inx P n (t), which yields 
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where 



f\x\ = -T D (e>*-l) 



(9) 



contains the counting field. Thus, the time-evolution of 
the CGF couples only to two other correlations (e mx J z ) 
and (e mx J 2 ). We therefore calculate the EoM for corre- 
lations of arbitrary order a > 1, which is given by 



1 <e^J«) = ]TTr {e^J«pW(t)} 



(10) 



r s n s [(J 2 + 2-l) Q - J«]X_ 



+r s (l + n s )[(J z -2-l) Q - J«]X+ 
+r D n Z5 [e-« ( J z + 2 • 1)" - Jf\ X 

+T D (l + n D ) [e ix (J z - 2 • 1)" - Jf] X H 



with 



fj 2 - J 2 ± 2 J, 



(11) 



and thus couples the evolution of (e mx J") to all powers 
of (e mx J^ where a' G {0, 1, . . . , a+2}. This leads to a 
hierarchy problem, such that the straightforward solution 
of these EoMs without factorizations is difficult for large 
j = Equations (8) - (11) are the first central result 
of this paper. 

For N = 1 (i.e., j = \) however, the above system can 
be solved exactly. For this special case we have J z = a z 
and (J 2 = 1 and thus we can solve the system of equations 
without the need of any factorization. We introduce the 



abbreviation A(x, t) = ( e \n X ) an d obtain 
±C(x,t) = -(?-*-!) ^-n D [A( X ,t)-l] 



_d 

dt 
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+ (e+ i *-l) 1 f[l + n D ] [A( X ,t) + l] 



izMx, t) = dt )- -A(x, t) ( i-Ax, t) ) , (12) 
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with 



-r s n s [A( Xl t) 
-— n D (e x 

~[l + n D ] (e +ix 



1]-T S [l + n s ] [A( X ,t) + l] 
1) [A( X ,t)-l] 

+ l)[A( X ,t) + l]. (13) 



The quadratic equation for A(x,t), cf. Eq. (12), can be 
solved and inserted in the equation for the CGF. Since 
we are interested in long-time dynamics, we can directly 
take derivatives with respect to x and send t — » oo to ob- 
tain the corresponding stationary cumulants. Note that 
to circumvent the quadratic equation one could alterna- 
tively calculate the dynamics of the moment-generating 
function A4(x,t) = (e mx ), where one has to solve a 
system of coupled linear equations and then construct 
the cumulants from the moments. As expected, the thus 
obtained CGF coincides in the long-time limit with the 
eigenvalue Ao(x) with smallest modulus, obtained by di- 
agonalizing the corresponding 2x2 Liouvillian, and is 
given by 



r s (l + 2ra< ? )+r D (l + 2V) (14) 



- \v s T D [e- ix (l + n s )n D + e +ix n s (l + n D )] 

(r s (l + 27i s )) 2 + (r D (l + 2n D ))] 
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In contrast, for large TV it is much more difficult to ob- 
tain the CGF, as already the size of the Liouvillian grows 
as (N + 1) x (N + 1). However, for the long-time limit of 
the first cumulant ((ii)) we can make use of Eq. (7) and 
the fact that trace conservation implies Tr {C(0)p} = 0, 
such that we find 

tME 



<</r>) = -iTrs{£'(0)p} 

= [ns - n D ) 

1 s 

(N — 2n) (1 + n) 



N+l 



cjn , (15) 
n N+i (2 + N + 2n) 



(1 + n) 



N+l 



n 



N+l 



where the superscript ME denotes the master equation 
solution. This result (which we obtained previously [16]) 
gives us one analytic benchmark for approximate solu- 
tions. 
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Due to [C'(x), £(x)) 0, the above approach yields no 
simple analytic results for higher cumulants (that require 
higher derivatives with respect to x)- The corresponding 
expressions still involve the full x-dependent Liouvillian, 
that needs to be diagonalized for a solution. To bypass 
the expensive numerical calculation, which becomes un- 
feasible already for moderate N, we have to apply addi- 
tional approximations. 



IV. APPROXIMATE EQUATIONS OF MOTION 

Eq. (8) and (10) constitute the EoM for the full CGF. 
To circumvent the problem of the hierarchy of operator 
equations, we introduce a lowest-order factorization at 
the level of the CGF and system operators, and show that 
already this very crude assumption gives useful results 
comparable to the master equation solution. 

The simplest way to make Eq. (8) solvable is to add 
one further assumption to the Born-Markov-Sccular ap- 
proximation, namely that the counting operator n has no 
correlations in any order a with the system operators 



(16) 



This assumption is strictly valid only when the statistics 
of particles counted at the drain and the system state 
are independent. Inserting the approximation Eq. (16) 
into Eq. (8) yields another central result of this work: an 
explicit equation for the CGF 



-C( X ,t)=n D (f[ x ]-f[-x])(J z ) 



fix] 



(f[x] + f[-x])[J 

[J 2 -(J 2 Z )+2{J Z )} , 



(17) 



(1,2) 



with f[x] gi vcn b Y E q- (9). 

Independent of the choice of solution for ( J- 
ing derivatives with respect to X shows that all 
derivatives of odd (((^2fe+i))) and even (((ri2fc))) cumu 
lants arc identical, respectively, and are given by 



tak- 
time 



d ((n 2 fc+i)) 

dt r D 

d ((n 2 fc)) 

dt 



n D 



[J 2 -{JI)+2{J., 
(J 2 - (J 2 ,)) 



>] , (18) 



2 (J z 



Therefore, provided that Eq. (16) is approximately valid, 
the statistics of the model can be retrieved with knowl- 
edge of only (J z ) and (J% V 



A. Approximate solutions for large spin 

For our specific system we are able to obtain exact 
steady state expressions for the expectation values of the 
system operator for arbitrary power (</"}, cf. Eq. (7). If 



one applied our approach to other systems, this would 
generally not be the case. To close the system of equa- 
tions wc would therefore have to calculate the EoM for 
the system operators as well, even when factoring corre- 
lations of the type of Eq. (10). 

To elucidate the validity of the central approximation 
Eq. (16), we also calculate solutions of Eq. (17) using 
approximate expressions for the evolution of the system 
operators. For this aim, using Eq. (3) yields an EoM for 
powers of J z given by 



dt r 



n([(J z + 2-l) a -J?]X_) 

+(n + l)([(J z -2.ir-J«]X + ) , (19) 



where T = T$ + Fjj, X± are given in Eq. (11) and 
n is given in Eq. (6). Since X± cx J^,J Z , the evolu- 
tion of an operator J" couples to all operators J z with 
a' G {l,2,...,a + l}. The system of coupled linear dif- 
ferential equations for (J z ) gives again rise to a hierarchy 
problem, now on the level of system operators, and can be 
solved approximately for large j. Factoring the expecta- 
tion value at any level (e.g. <^ Q+ "') « (J?) W*'}) wiU 
lead to a coupled set of non-linear differential equation, 
that in the time-dependent case requires numerical solu- 
tions. In our approach we are interested in steady-state 
results and will therefore always use long-time expecta- 
tion values. 



Approximation 1 

To estimate the errors introduced by the factorizations 
on top of the general one (16), we can use solutions of 
Eq. (17} as a benchmark, where the exact stationary re- 
sults (J z ) ex an d (Jf) °f Eq. (7) for a = 1, 2 are in- 
serted. Since we have applied our general factorization 
Eq. (16) even when using the exact expressions we will 
denote this solution as Approximation 1. In the following 
we differentiate between two further levels of factoriza- 
tion. Higher order approximations can easily be con- 
structed. 



Approximation 2 

Here, in addition to (16) we take only the dynamics of 
(J z ) into account, assume (Jf ) ~ (J z ) in Eq. (19) with 
a = 1 (Approximation 2). This is valid whenever the 
variance vanishes and leads to a single quadratic equa- 
tion, 



^^--(^+l)(Jz) + l(Jz) 2 



1 



(20) 



which can be solved analytically for all times. Since 
we are interested in steady state dynamics of the FCS, 
we take the long-time limit and insert the solution in 
Eq. (17) where we again assume (J z ) ~ (Jz) 2 - 
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Approximation 3 

Alternatively, again in addition to (16) we take (</ 2 ) hi 
Eq. (19) for a = 1 into account and factorizc third order 
correlations (J 3 ) w (J z ) for Eq. (19) with a = 2 
(Approximation 3). This yields a coupled system of (non- 
linear) differential equations, 



d (J z ) 

dt r 
dt r 



(2n + l)(J„) 



2 x 7 2 



(21) 



= (2n+l)[j 2 -3(j 2 )] - [J 2 -2] (J,) 



that in general has to be solved numerically. However, 
the steady state solution can be obtained analytically and 
again is inserted into Eq. (17), where apart from Eq. (16) 
we now do not have to assume further factorizations. 



B. Thermodynamic limit for the 
cumulant-generating function 

The full CGF for the master equation solution is not 
accessible for arbitrary system size, since the diagonal- 
ization of the large Liouvillian is computationally cum- 
bersome. For the approximate case however, solutions 
can be easily derived by solving Eq. (17) with the differ- 
ent approximations and we are thus able to examine the 
thermodynamic limit of infinite system size. 

The obtained expressions for the approximate CGFs, 
with a marking the applied approximation, arc 
C a ( X ,t) = F( X ,t)C a with 



C 1 



(fi + l)(JV-2n)+n^J (A + 2n + 2) 



(l+n) 



N 



C 2 = - (2n + 1) + v/(2n + l) 2 + JV(2 + N) , 
C 3 = -i (3 + 6n - (2n + 1) _1 ) 

+J\ (3 + 6n - {2n + l)" 1 ) 2 + A^(2 + N) , (22) 
where the common function 



F{ X ,t)=T D t 



He 



(e ix - l)(n D + l)n 



l)n D (n + l) 



(23) 



is shared by all approximations. The form of Eq. (23) 
shows that the transport process is a balance of drain oc- 
cupation and effective average occupation, and the choice 
of approximation yields an overall scaling of the corre- 
sponding emission and absorption processes. 

These expressions allow us to directly verify the ther- 
modynamic limit for N — > oo in specific limits. 



For n <C N we have a linear scaling of the CGF in the 
system size and obtain the same result for all approxi- 
mations 



1 



N 



lim C a ( X ,t) = F( X ,t). 



N 



(24) 



If we perform the thermodynamic limit and simulta- 
neously require n 3> N, we are always in the super- 
transmittance regime [16]. Scaling only the source oc- 
cupation, with a finite drain occupation, we recover the 
quadratic scaling of the CGF on the system size for all 
approximations 



lim C a ( X ,t) 

ns—^oo 



-N(N + 2) 



1) n D 



t. 



l)(l + n D ) 



(25) 



where the prefactor gives a deviating scaling for Approx- 
imation 2 (x2 = 4) but agrees with the exact solution 
(xi = x% = 6) for the high thermo-bias, cf. [16], other- 
wise. This indicates that the coherent effect of "super- 
transmittance" [16] is correctly found also by the factor- 
ization approach. 

As a further benchmark for our approximations we 
check the limit of low thermo-bias, where for low source 
occupations rig <S 1 and vanishing drain occupations 
Ud = we again exactly recover previous findings 
(cf. [16]) for the master equation solution 



c(x,t) 



S 1 D 



with all approximations. 



(e +ix - l)n s Nt, 



(26) 



C. Scaling of stationary cumulants 

To quantify the errors which are introduced by the 
different approximations, we compare the corresponding 
long-time cumulants to the full master equation solution, 
which were derived in [16]. 

First, we compare the ratio of approximate and master 
equation solution for the steady state of the first current 
((i AP )) 

cumulant ^jfceyy- The master equation solution is given 

in Eq. (15) and the approximate solutions are calculated 
with Eq. (22). 

Interestingly, approximation 1 yields the exact analytic 
expression for arbitrary source or drain occupations n a , 
tunneling rates T a and system size N for the first sta- 
tionary current cumulant, as shown by the black, solid 
line in all panels of Fig. 2. 

For low source occupation ns, cf. top panel of Fig. 2, 
approximations 2 (red, dotted line) and 3 (blue, dashed 
line) yield results comparable to the exact solution. For 
higher source occupation ns (middle and bottom panel), 
approximation 2 shows a deviation from the exact results 
for small system size. With higher occupations it is sus- 
tained up to larger system sizes. Approximation 3 shows 
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FIG. 2. (Color Online) Logarithmic plot of the ratio of ap- 
proximated and exact master equation solution for the first 

steady state cumulant ^ 7 me^ versus the system size N for 

the different approximations 1 (grey filled, black circles), 2 
(red squares), 3 (filled, blue diamonds), cf. Eq. (22), com- 
puted for different dimensionless source occupations ns =0.1 
(top), ns = 1 (middle) and ns = 10 (bottom). Note the dif- 
ferent scales on the ordinate axes. Other parameters: no — 0. 



FIG. 3. (Color Online) Double logarithmic plot of the second 
steady state cumulant {{I2)) in units of T = Fg = To ver- 
sus the dimensionless source occupation ns for the numerical 
master equation solution (bold dot-dashed, black) and the 
different approximations 1 (solid, black), 2 (dotted, red), and 
3 (dashed, blue), cf. Eq. (22), computed for different system 
size N — 5, 10, 20, 40, 80 (bottom to top). Other parame- 
ters: n_o = 0. 



a different type of deviation. The position of the maxi- 
mal deviation is dependent on the source occupation and 
shifted to larger system size for higher occupations. Ad- 
ditionally, as expected all approximations approach the 
exact solution with larger system sizes. 

The approximate solutions for the second cumulant 
((i^ p )} can be analytically calculated from Eq. (22). 
However, for the master equation results ((l2 IE )), wc 
have to revert to numerical solutions. Since we are in- 
terested in the validity of our factorization Eq. (16) for a 
broader range of occupations, we compare in Fig. 3 the 
approximated second cumulants ((l£ p )) and the numeri- 
cal master equation solution ((l2 IE )) for different system 
sizes N = {5, 10, 20, 40, 80} over the range of source oc- 
cupations ns 6 {0.01 . . . 1000} (setting no = 0). 

While for low source occupations all approximations 
show perfect agreement with the master equation solu- 
tion (thick, black dot-dashed line), approximation 2 (red, 
dotted line) shows a constant deviation to the master 
equation solution for large source occupations. However, 
approximations 1 (thin, black solid line) and 3 (blue, 
dashed line) show excellent agreement, again. This is 
valid for all considered system sizes. For intermediate 
source occupations all approximations show large devia- 
tions from the master equation solution in higher cumu- 
lants. 



D. Quality of approximations 

The quality of the approximations can be under- 
stood by evaluating the steady state expression (J"), cf. 



Eq.(7), for a = 1, 2, 3 in the corresponding regimes for 
n — > or n — > 00. In the limit of vanishing source occu- 
pation n s -> 0, both ( Jf) « {J z f and ( J 2 3 ) « (J 2 Z ) (J z ) 
become exact. In the limit of infinite source occupation 
ns — > 00 however, all odd expectation values {J z a+1 ) 
are vanishing, while even expectation values are scaling 
with powers of the system size N. Thus, (J z ) ~ (J z ) 2 is 
not a good approximation in this regime, while (J z ) ~ 
(J z ) (J z ) again gets exact, explaining the failure of ap- 
proximation 2. 

The validity of the general approximation Eq. (16) can 
be clarified by calculating 

N 

OO — 

(e^J«)= E (n,m\e inx J:p\n,m) (27) 

tl=-M m =-| 
N 

CO "2" 



J2 E e™*(m|J>W 



n=-oo m —_ 



J2 E ( 2m T einx H p [n) \ m ) . 



n=-co m= _N 

where = (n\p\n) is the conditioned density ma- 
trix, and in the last line the eigenvalue of J" was in- 
serted. Obviously, the expectation value can be fac- 
tored when the system is close to the ground state 
(m| \m) ~ S m _«.C n , which is the case in the limit 

n — > (e.g. low thermobias and no = 0). Alternatively, 
factorization is possible when the system is equiparti- 
tioned (m \ p^ |ra) ~ C n , which is fulfilled for the limit 
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n — >• oo (e.g. large thcrmobias). This explains the good 
performance of the approximated solutions in the regimes 
of low and high weighted average occupations n. 

In contrast, in the intermediate region of n neither 
fluctuations are suppressed due to the system being close 
to its ground state, nor the supply of bosons from the 
source is large enough to re-pump the large spin J z to 
its thermalized state. Thus, the dynamics is governed by 
higher order correlations in the intermediate regime that 
are not taken into account at this level of factorization. 



V. CONCLUSION 

The application of the equation of motion method to 
the CGF on top of the master equation allows for the 
approximate calculation of the full dynamics of the sys- 
tem. For the special case N = 1 the method becomes 
exact and recovers the FCS results obtained by diagonal- 
ization of the full Liouvillian. In the limits of low and 
high thcrmo-bias the master equation results of previous 
work [16] are recovered for the full long-term CGF. Fur- 
thermore, this method enables us to directly access the 
thermodynamic limit, which is not generally possible in 
the master equation solution. 

For the approximate solutions we showed that it is nec- 



essary to take at least the equation for (J%\ into account 
and factorize at the level of to get comparable re- 

sults to the master equation in the stationary limit. 

Future work should be directed to the question whether 
the inclusion of higher correlations of (e mx J") into the 
equation for the CGF yields better results for higher cu- 
mulants and also at intermediate source/drain occupa- 
tions. 

As is evident from the case of N = 1, one should in 
general carefully examine whether for system operators 
A our main approximation (e mx 'A) w (e mx ) (^4) is suf- 
ficient, or if one has to close the system of equations 
at higher order. However, since the zeroth order of the 
method allows for closure with only a few equations, and 
many relevant probability distributions are close to Gaus- 
sian and thus governed by the first two cumulants [27], 
this approach could find application in a broader com- 
munity. 

ACKNOWLEDGMENTS 

We acknowledge support by the DFG via GRK 1558 
(M.V.), grants SCHA 1646/2-1 (G.S.), BRA 1528/7, 
BRA 1528/8, SFB 910 (T.B.). We have benefited from 
discussions with D. Braun. 



[1] Y. M. Blanter and M. Biittiker, Phys. Rep. 336, 1 (2000). 
[2] G. Ferrini, A. Minguzzi and F. W. J. Hekking, Phys. Rev. 

A 80, 043628 (2009). 
[3] J. Gabelli andB. Reulet, Phys. Rev. B 80, 161203 (2009). 
[4] R. J. Cook, Phys. Rev. A 23, 1243 (1981). 
[5] F. Schlogl and E. Scholl, Z. Phys. B 51, 61 (1983). 
[6] L. S. Levitov, H. Lee and G. B. Lesovik, J. Math. Phys. 

37, 4845 (1996). 
[7] D. A. Bagrets and Y. V. Nazarov, Phys. Rev. B 67, 

085316 (2003). 

[8] C. Flindt et al, Proc. Natl. Acad. Sci. U.S.A. 106, 10116 
(2009). 

[9] A. Andreev and A. Kamenev, Phys. Rev. Lett. 85, 1294 
(2000). 

[10] M. Richter, A. Carmele, A. Sitek and A. Knorr, Phys. 

Rev. Lett. 103, 087407 (2009). 
[11] R. H. Dicke, Phys. Rev. 93, 99 (1954). 
[12] M. Gross and S. Haroche, Phys. Rep. 93, 301 (1982). 
[13] R. Bonifacio, P. Schwendimann, and F. Haake, Phys. 

Rev. A 4, 302 (1971). 
[14] F. Haake and R. J. Glauber, Phys. Rev. A 5, 1457 (1972). 
[15] P. A. Braun, D. Braun, F. Haake, and J. Weber, Eur. 

Phys. J. D 2, 165 (1998). 
[16] M. Vogl, G. Schaller, and T. Brandes, Ann. Phys. (NY.) 

326, 2827 (2011). 



[17] T. Brandes, Phys. Rep. 408, 315 (2005). 
[18] R. H. Hadfield, Nat. Phot. 3, 696 (2009). 
[19] S. De Liberato, N. Lambert, and F. Nori, Phys. Rev. A 

83, 033809 (2011). 
[20] I. S. Grudinin, H. Lee, O. Painter, K. J. Vahala, Phys. 

Rev. Lett. 104, 083901 (2010). 
[21] J. G. Bohnet et al, Nature 484, 78 (2012). 
[22] G. Schaller, G. Kiefilich and T. Brandes, Phys. Rev. B 

80, 245107 (2009). 
[23] H.-P. Breuer and F. Petruccione, The Theory of Open 

Quantum Systems, (Oxford University Press, Oxford, 

2002). 

[24] G. Schaller, Phys. Rev. E 83, 031111 (2011). 

[25] G. S. Agarwal, Quantum Statistical Theories of Sponta- 
neous Emission and their relation to other approaches, 
Springer Tracts in Modern Physics Vol. 70 (Springer, 
Berlin, 1974). 

[26] C. Flindt, T. Novotny, A. Braggio, M. Sassetti and A.P. 

Jauho, Phys. Rev. Lett. 100, 150601 (2008). 
[27] T. Karzig and F. von Oppen, Phys. Rev. B 81 045317 

(2010). 



